Identification of an immune‐related gene signature as a prognostic target and the immune microenvironment for adrenocortical carcinoma

Abstract Background Adrenocortical carcinoma (ACC) is a rare endocrine malignancy. Even with complete tumor resection and adjuvant therapies, the prognosis of patients with ACC remains unsatisfactory. In the microtumor environment, the impact of a disordered immune system and abnormal immune responses is enormous. To improve treatment, novel prognostic predictors and treatment targets for ACC need to be identified. Hence, credible prognostic biomarkers of immune‐associated genes (IRGs) should be explored and developed. Material and methods We downloaded RNA‐sequencing data and clinical data from The Cancer Genome Atlas (TCGA) data set, Genotype‐Tissue Expression data set, and Gene Expression Omnibus data set. Gene set enrichment analysis (GSEA) was applied to reveal the potential functions of differentially expressed genes. Results GSEA indicated an association between ACC and immune‐related functions. We obtained 332 IRGs and constructed a prognostic signature on the strength of 3 IRGs (INHBA, HELLS, and HDAC4) in the training cohort. The high‐risk group had significantly poorer overall survival than the low‐risk group (p < .001). Multivariate Cox regression was performed with the signature as an independent prognostic indicator for ACC. The testing cohort and the entire TCGA ACC cohort were utilized to validate these findings. Moreover, external validation was conducted in the GSE10927 and GSE19750 cohorts. The tumor‐infiltrating immune cells analysis indicated that the quantity of T cells, natural killer cells, macrophage cells, myeloid dendritic cells, and mast cells in the immune microenvironment differed between the low‐risk and high‐risk groups. Conclusion Our three‐IRG prognostic signature and the three IRGs can be used as prognostic indicators and potential immunotherapeutic targets for ACC. Inhibitors of the three novel IRGs might activate immune cells and play a synergistic role in combination therapy with immunotherapy for ACC in the future.


| INTRODUCTION
Adrenocortical carcinoma (ACC) is a rare endocrine malignancy worldwide. 1,2 The incidence is reported to be one to two per million per year. The prognosis is poor with a median overall survival (OS) of 3-4 years. 3 Even with complete tumor resection and adjuvant therapies, the prognosis of patients with ACC following treatment remains unsatisfactory. 4 The risk of recurrence and metastasis ranged from at least 15% to 40% among patients with ACC after initial treatment. 5,6 Consequently, novel prognostic predictors and treatment targets for ACC are needed and credible prognostic biomarkers of immune-associated genes (IRGs) should be explored.
In the microtumor environment, the impact of a disordered immune system and abnormal immune responses is enormous. 7 Immune cells may modify the signaling molecule to influence the immunological response, which significantly impacts the risk of recurrence, change, and diffusion of the neoplasm. 8,9 Previous research studies have reported that some tumor cells may evade immune surveillance and immune elimination, resulting in tumor infiltration and metastasis. 10 An abnormal immune state is thought to be associated with glioblastoma progression. 11 Martinez-Bosch et al. 10 suggested that immunosuppressive macrophages and myeloid-derived suppressive cells may significantly suppress the immune microenvironment in pancreatic cancer. Peng et al. 12 identified two immune-related genes could serve as potential biomarkers for immunotherapy in ACC, which provide better insights into ACC microenvironment. In addition, N6-methyladenosine methylation regulators had impacts on ACC prognosis through regulating immune-related functions. 13 Bioinformatics analyses based on The Cancer Genome Atlas (TCGA) constructed a competitive endogenous RNA network associated with tumor-infiltrating immune cells and determined that the immune system is implicated in the microenvironment of ACC. 14 Aberrantly expressed IRGs may be potential therapeutic targets and prognostic biomarkers for patients with ACC. For this purpose, we explored IRGs and constructed an IRG prognostic signature. This study was designed using the gene expression profiles from the TCGA ACC and Genotype-Tissue Expression (GTEx) data sets. To investigate the IRG signature's feasibility, we validated the predictive power in two Gene Expression Omnibus (GEO) data sets. We also evaluated the correlation between the IRG signature and tumor-infiltrating immune cells.

| TCGA ACC and GTEx data sets
The TCGA ACC data set provided ACC specimens (n = 79), whereas the GTEx data set provided normal adrenal specimens (n = 258). We acquired the RNAsequencing (RNA-seq) data and clinical samples from the GTEx (http://www.gtexportal.org/home/index.html) and TCGA (http://portal.gdc.cancer.gov/) databases. Differentially expressed genes (DEGs) were aberrated transcripts identified under different biological context through RNA-seq analysis. 15 DEGs makes RNA-seq analysis possible to get a deep understanding of pathogenesis-related molecular mechanisms and biological functions. Therefore, differential analysis has been regarded as valuable for diagnostics, prognostics, and therapeutics of tumors. To identify DEGs from the ACC and normal adrenal samples, we used the edgeR R package. 16 The cutoff criteria were |log Fold Change| ≥ 1 and an adjusted p (adj. p) < .05.

| Identification of an immune-related gene prognostic signature
We construct the IRG signature by using the survival R package. The identification and verification criteria for IRG signatures in the ACC specimens were as follows: (1) intact information regarding the IRG expression and clinical features (age, gender, tumor TNM stage, and survival time); and (2) specimens with total survival <30 days were excluded (some nonneoplastic factors such as severe infection or hemorrhage may have destroyed these specimens). From these, we selected 77 ACC samples for the IRG signature construction. The 77 specimens were randomly divided into the training (n = 39) and the testing (n = 38) cohorts. The training cohort was utilized for the IRG signature construction. The other cohort was utilized for validation of the IRG signature.
Next, we screened for prognosis-related IRGs (p < .01), regardless of whether the IRGs were protective or deleterious according to hazard ratio (HR) as well as 95% confidence interval (CI). The 332 IRGs were analyzed by univariate Cox regression analysis in the training cohort. After that, we selected an optimal model from the prognosis-related IRGs using the Akaike information criteria (AIC) method. The three IRGs with minimum AIC values were finally selected for the prediction model construction.
We analyzed the prognosis-related IRGs using multivariate Cox regression analyses to construct a prognostic IRG signature and identify the coefficients. [18][19][20] The following in-house formula was used to calculate the risk score of the prognostic IGR signal for each specimen: Risk score = Expression IRG1 × Coefficient IRG1 + Expression IRG2 × Coefficient IRG2 + … + Expression IRGn × Coefficient IRGn . The risk score of the prognostic IRG signature may be determined based on a linear combination of the IRG expression level weighted by the regression coefficients. We acquired the coefficient by log-transformed HR, derived from the multivariate Cox regression analysis. 21,22 We divided the specimens into the low-risk and high-risk groups based on the median value of the risk score. Univariate and multivariate Cox regression analyses, time-dependent receiver operating characteristic (ROC) curve, and Kaplan-Meier (KM) survival curve were conducted to evaluate the predictive value of the risk score in the training cohort.

| Internal validation and external validation of the IRG signature
To further validate the predictive ability of the IRG signature, we used the testing cohort and the entire F I G U R E 1 Workflow of this study. The study was carried out in The Cancer Genome Atlas (TCGA) adrenocortical carcinoma (ACC), Genotype-Tissue Expression (GTEx), and Gene Expression Omnibus (GEO) data sets. Differentially expressed genes (DEGs) were calculated between ACC samples from TCGA ACC data set and normal adrenal samples from GTEx data set. GSEA was conducted based on all DEGs and we found the DEGs were enriched in immune-related functions. Then, 322 immune-associated genes (IRGs) were extracted from Molecular Signatures Database v4.0. The training cohort was used to identify prognostic IRGs and establish a prognostic signature based on three IRGs (INHBA, HELLS, and HDAC4). The prognosis analysis was validated in the testing cohort and the entire TCGA ACC cohort, respectively. In addition, external validation was carried out based on GSE10927 cohort and GSE19750 cohort. Tumor-infiltrating immune cells analysis was performed based on CIBERSORT tool to investigate the association between the three-IRG signature and immune system. TCGA ACC cohort to perform internal validation, and used GSE10927 and GSE19750 cohorts from GEO to perform external validation. The validation was conducted by following the same analysis methods of the training cohort.

| Functional enrichment analysis between different risk groups and correlation analysis
Principal component analysis (PCA) was carried out to profile the expression patterns of the low-risk and high-risk groups according to the IRG signature by using scatterplot3d R package. We identified DEGs in both groups and adopted the clusterProfiler and GOplot R packages for GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses 23 to uncover the potential capacities in the two risk groups and the TCGA ACC cohort. Furthermore, the adj. p < .05 was used as the cutoff value. We investigated the pairwise gene correlation of three IRGs (   We retrieved 22 types of tumor-infiltrating immune cells from CIBERSORT, a web tool for analyzing tumor-infil trating immune cells (http://cibersort.stanford.edu/). 24,25 To identify the differences between the low-risk and high-risk groups, we measured the quantity of each immune cell. p < .05 was set as the cutoff value.

| Statistical analysis
Statistical analysis was carried out by GraphPad Prism 7.0 and R software (v3.5.3: http://www.r-project.org). The RNA-seq data were log2-transformed. To assess the influences of the IRG prognostic signature on OS and other clinical features, Cox regression analyses, log-rank test, and KM method were used. The χ 2 test and Student's t test were adopted to measure qualitative variables and quantitative variables, respectively. According to the IRG signature and other indexes, we evaluated the predictive value using time-dependent ROC. The cutoff value of the two-sided p was set at .05 (two-sided p < .05). Figure 1 shows the study flowchart.

| Identification of DEGs and functional enrichment analysis
All genes from the TCGA ACC and GTEx data sets were analyzed via gene differential expression analysis. From these, we selected 79 ACC specimens and 258 normal adrenal specimens. We used edgeR R package to identify DEGs between ACC and normal adrenal samples to further explore the immune-related DEGs and immunerelated functions. In total, we identified 5609 DEGs from these specimens, including 3519 upregulated and 2090 downregulated DEGs (Figure 2A). The biological capacities of DEGs determined via GO enrichment analysis was identified by GSEA. The results indicated that upregulated genes were involved in the humoral immune response (GO:0006959; p = .027; Figure 2B). Moreover, downregulated genes correlated with the effector process of the immune system (GO:0002252; p = .026) and immune response (GO:0006955; p = .038; Figure 2C,D), which showed that the tumorigenesis and progress of ACC are related to the immune system and immune responses.

| Identification of immune-related genes and the clinical characteristics of TCGA patients with ACC
We identified and obtained 332 IRGs from the Molecular Signatures Database v4.0. To construct and validate the IRG signature, 77 patients with ACC were included. These patients were randomly divided into the training (n = 39) and the testing (n = 38) cohorts. Table 1 was the clinical features of 77 ACC patients.

| Identification of the IRG prognostic signature in the training cohort
We identified prognostic IRGs from the 262 IRGs in the training cohort using univariate Cox regression analysis. There were seven prognosis-related IRGs.
To construct an IRG signature, a multivariate Cox regression model with the smallest AIC was applied to the seven prognostic IRGs. Ultimately, three IRGs were selected and we constructed the three-IRG prognostic signature. 3.4 | Assessing the predictive capacity of the three-IRG prognostic signature for patients with ACC in the training cohort From the training cohort, 39 patients with ACC patients were divided into two groups based on the median value of the three-IRG signature risk score, namely, the high-risk (n = 19) and low-risk (n = 20) groups. Figure 3A,B summarizes the details of each patient in the training cohort, including survival time, survival status, and the risk score. Figure 3C T A B L E 2 The three immune-related genes identified from Cox regression analysis  illustrates the expression of the three IRGs. The KM survival curve of all the patients with ACC in the training cohort is shown in Figure 3D. We analyzed the survival time and found that patients in the lowrisk group had a longer survival time than those in the high-risk group (p < .001). As shown in Figure 3E, by plotting time-dependent ROC curves of the three-IRG prognostic signature and other clinical features in the training cohort, we identified that the area under the ROC curve (AUC) value was 0.923. We also evaluated the age, three-IRG prognostic signature, gender, and tumor stage using multivariate Cox regression analyses, as shown in Figure 3F. From the risk score (HR = 2.585, 95% CI = 1.569-4.258, p < .001) of the three-IRG signature, a unique prognostic target of overall poor survival was identified (Table 3).
3.5 | Authentication of the three-IRG prognostic signature in the entire TCGA ACC cohort and the testing cohort We utilize the three-IRG signature to evaluate the OS of patients with ACC in the testing cohort to prove the stability and predictive capacity of the three-IRG signature. Risk scores of the three-IRG signature in 38 patients with ACC from the testing cohort were computed using the previously mentioned formula. The 38 ACC patients were divided into the high-risk and low-risk groups (n = 19 per group) based on the median value. Figure 4A-B provides the clinical features of the testing cohort, the survival time, survival status, and the risk score. The expression of the three IRGs is shown in Figure 4C. Figure 4D displays the KM survival curve of the testing cohort. Our results found that patients of the lowrisk group subsisted for a longer time compared with the high-risk group (p = .002). Figure 4E shows the plotting time-dependent ROC curves of the three-IRG prognostic signature and other clinical features in the training cohort, which demonstrated that the AUC of the three-IRG prognostic signature was distinctly better than other clinical features (risk value was 0.867). The results between the training cohort and the testing cohort were the same. From the multivariate Cox regression analysis, the risk score (HR = 2.168, 95% CI = 1.281-3.668, p = .004) confirmed the three-IRG signature as a unique prognostic target ( Figure 4F and Table 3).
All the patients in the TCGA ACC cohort were analyzed using time-dependent ROC curves, KM survival curve, and multivariate Cox regression analysis. The results are shown in Figure 5 and Table 3.

| External validation using the GSE10927 and GSE19750 cohorts
To externally validate the three-IRG signature in the GSE10927 (n = 24) and GSE19750 (n = 22) cohorts, we applied the time-dependent ROC curve (Table 4), univariate and multivariate Cox regression analyses, and KM survival curve ( Table 3). The KM survival curve and time-dependent ROC curve for the GSE10927 cohort are shown in Figure 6A,B. The results confirmed that patients of the low-risk group subsisted for a longer time compared to those of the high-risk group (p = .005). The AUC of three-IRG prognostic signature was distinctly better than other clinical features (risk value was 0.861). Similar results were found in the GSE19750 cohort ( Figure 6C,D).

| Survival, clinical features, and the three-IRG prognostic signature
We analyzed the entire TCGA ACC cohort by stratified survival analysis, to further verify the prognostic capacity and explore the extensive feasibility of the three-IRG signature ( Figure 7A-G). The samples with Stage I-II (p = .011), Stage III (p = .008), and Stage IV (p = .031) were assigned to the low-risk group, which had better survival than the high-risk group. Similar results were found in diverse ages and different gender. Clinical features were distributed along with the three-IRG prognostic signature. The high-risk group had a higher number of mortalities, which was reflective of poorer survival in these patients who had a high-risk value (p < .05) ( Figure 8B-D). Furthermore, there were different risk scores for patients in different disease stages. Patients at the inchoate stages (Stage I and Stage II) have lower risk values than those at more advanced stages (Stage III and Stage IV; p < .05; Figure 8E-G), which verified that risk scores of the IRG prognostic signature were distinctly related to the development of ACC.

| PCA and functional enrichment analysis
To determine the disparate distribution patterns between the high-risk and low-risk groups based on the three IRGs, we performed a PCA analysis. The two groups were obviously separated in two different directions based on the three IRGs. It was evident that the samples in the high-risk group were significantly different from that of the low-risk group ( Figure 9A-C). From the entire TCGA ACC cohort, 664 upregulated and 563 downregulated DEGs were revealed between low-risk and high-risk groups. We identified their capacity using the KEGG and GO functional enrichment analyses. The top 15 terms of GO are listed in Figure 9D and Table 5. The top 15 terms of KEGG are shown in Figure 9E and Table 6. We found that they were significantly enriched in the extracellular matrix regulation, cell cycle, and phosphatidylinositol 3-kinase (PI3K)-Akt signaling pathway.

| Survival analysis and correlation analysis of INHBA, HELLS, and HDAC4
We further explored the prognostic value of the three IRGs (INHBA, HELLS, and HDAC4) for construction of the signature and identified that higher expression levels of INHBA, HELLS, and HDAC4 were related to worse OS time among patients with ACC (p < .05; Figure 10A-C). Higher expression levels of HELLS and HDAC4 also indicated worse disease-free survival ( Figure 10D-F). INHBA, HELLS, and HDAC4 in the TCGA ACC data set were analyzed by pairwise correlation analysis. The expression levels of HELLS and HDAC4 demonstrated an obvious positive correlation (p < .001 and R = 0.648). An increase in HELLS was related to an increase in HDAC4 ( Figure 10G).

| Correlation between tumorinfiltrating immune cells and the three-IRG signature
We compared the tumor-infiltrating immune cells between low-risk and high-risk groups in the entire TCGA ACC cohort, to investigate the relationship between tumor immune microenvironment and the three-IRG prognostic signature. T cells, natural killer (NK) cells, macrophage cells, mast cells, and myeloid dendritic cells showed disparate abundance in the lowrisk and high-risk groups (p < .05; Figure 11). The correlation heat map was also plotted based on the distribution of immune cells in patients with ACC ( Figure 12).

| DISCUSSION
Accumulating evidence indicated that aberrantly expressed genes in tumors may be potential biomarkers for diagnosis, target therapy, and prognosis through bioinformatics analysis. 26 In cancer treatment, immunotherapy is gaining more and more attention. Currently, we can identify potential prognostic targets, analyze underlying mechanisms, and explore valuable IRGs through high-throughput sequencing data. 27,28 We researched IRGs using the TCGA ACC data set and constructed a three-IRG prognostic signature that may be used as a prognostic indicator or immune therapeutic target in patients with ACC.
The immune system impacts the tumor microenvironment, particularly tumor immune escape. 26 Immune components have quite a significant influence on the clinical outcomes and gene expression by tumor T A B L E 5 GO functional enrichment analysis between high-risk and low-risk groups in entire TCGA ACC cohort tissues, which exist in the tumor microenvironment. 9,29,30 There were 3519 upregulated and 2090 downregulated DEGs found in ACC samples and normal adrenal samples in this study. GSEA suggested that these DEGs are involved in the humoral immune response (GO:0006959), and the immune effector process (GO:0002252) and immune response (GO:0006955), which indicated that the development of ACC is significantly correlated with immune responses and the immune system ( Figure 2). We extracted 332 IRGs and constructed an IRG signature based on three IRGs (INHBA, HELLS, and HDAC4) in the training cohort. Based on the median value of the three-IRG signature risk score, we divided the samples into the low-risk and high-risk groups. We discovered that the prognosis of samples could be distinguished by the three-IRG signature. The OS of the high-risk group was worse than those of the low-risk group ( Figure 3D). The multivariate Cox regression analysis suggested that the risk score may be regarded as a unique prognostic target of OS through multivariate Cox regression analyses and the three-IRG signature was superior to other clinicopathologic features according to the ROC curve. The above findings were verified the testing cohort, and the whole TCGA ACC, GSE10927, and GSE19750 cohorts.
We then explored the wide applicability of the three-IRG signature. Survival analyses were conducted in disparate subgroups, stratified by age, tumor stage, and gender. We found that the three-IRG signature was applicable in various subgroups ( Figure 7). Moreover, we observed that the risk score of patients with early-stage ACC (Stage II and Stage I) was lower than the advanced one (Stage IV and Stage III; Figure 8). Our findings demonstrated that the three-IRG prognostic signature is not only a potential predictive indicator of tumor progression but can be applied to all patients with ACC.
Functional enrichment analyses demonstrated that the DEGs between low-risk and high-risk groups were obviously enriched in cell cycle function and PI3K-Akt signaling pathway. Cell cycle regulation is known to play a vital role in the differentiation and growth of tumor cells. 31 Subramanian et al. 32 identified a series of genes involved in regulating pathways of the cell cycle and overexpression of these genes was correlated with poorer OS in patients with ACC. LncRNA-HOTAIR participated in ACC development by regulating the cell cycle and proliferating ACC cells. 33 PI3K-Akt signaling pathway is a key driver in carcinogenesis and Akt overactivation has been verified in various endocrine gland neoplasms, including thyroid carcinoma subtypes, parathyroid carcinoma, pituitary tumor, and pheochromocytoma. 34 In vitro experiments revealed that low expression of SCTR could stimulate proliferation via the PI3K/AKT signaling cascade in ACC cells. 35 Recently, with the advent of checkpoint inhibitors and targeted immunotherapy, a greater focus has been placed on tumor immunotherapy for patients with ACC. 36,37 Our study evaluated the prognostic value of three IRGs for the construction of the IRG signature. Higher expressions of INHBA, HELLS, and HDAC4 were correlated with worse OS time for patients with ACC. Hence, the three IRGs may be promising therapeutic targets for ACC.
INHBA encodes a member of the transforming growth factor-β super-family of proteins. 38 INHBA participates in the regulation of the immune system, especially in skin suffering from repetitive ultraviolet-B irradiation. 39 Mesenchymal stromal cells inhibited monocytes through upregulating INHBA in patients with myelodysplastic syndrome. 40 IDO inhibitor induced the upregulation of IHNBA in regulating the functions of NK cells among sarcoma patients. 41 HELLS belongs to the SNF2-family, such as chromatin remodeling proteins. 42 HELLS mutation leads to immune deficiency syndrome in children. 43 A study based on adult mice suggested that although HELLS is not indispensable for lymphoid development, it is still necessary for the proliferation of peripheral T lymphocytes. 44 In addition, high-dose-rate γ-irradiation could impair immune signaling pathways by inducing the downregulation of HELLS. 45 HDAC4 is one of the histone-modifying enzymes, which are the main epigenetic regulators in the control of inflammatory processes. Furthermore, overexpression of HDAC4 may suppress the production of Type I interferons induced by pattern-recognition receptors in innate immunity. 46 Aberrant epithelial responses may arise via a self-defense mechanism generated by miR-22 through suppressing HDAC4 expression in epithelial cells. 47 Moreover, the HDAC4/PGRN and HDAC4/ nuclear factor-κB axes have been involved in the regulation of inflammatory cytokines in rheumatoid arthritis. 48 To identify gene signatures and further create clinical predictive models, high-throughput sequencing data has been utilized in large-scale research studies. Xu et al. stated that a signature could forecast the prognosis of ACC according to survival-associated alternative splicing events. 49 We constructed an immune-related signature based on 3 IRGs in the TCGA ACC data set and further validated its application in two GEO data sets. Furthermore, the predictive capacity of the signature was further confirmed. Besides, we also explored the association between tumor-infiltrating immune cells and the three-IRG signature and we found that T cells, NK cells, macrophage cells, myeloid dendritic cells and mast cells showed diverse abundance between the lowrisk and high-risk groups, which further verified the immune correlation of this signature. 14 Although the three IRGs were regarded as therapeutic targets for the treatment of ACC, the underlining mechanisms have not been clarified, which was one limitation of this study. Therefore, more in vitro and in vivo experiments based on a larger sample size are necessary for validating these findings and applying our results into clinical immunotherapy practice in the future.

| CONCLUSION
We identified an immune-related signature based on 3 IRGs in ACC. The three-IRG signature has important clinical significance, not only as a good classifier in distinct subgroups of ACC, but also as a unique prognostic target for ACC. Inhibitors of the three novel IRGs (INHBA, HELLS, and HDAC4) might activate immune cells in ACC immune microenvironment and play a synergistic role in combination therapy with immunotherapy for ACC in the future.